Retour à la page d’accueil


Lors de cette séance nous allons apprendre:

  1. A ajouter des résultats de tests statistiques à des graphiques

  2. A ajouter des statistiques descriptives à des graphiques

  3. A selectionner des lignes avec la fonction filter

  4. A représenter des données temporelles (gather et spread)


Préparation de l’environnement de travail

Tout d’abord ouvrez le projet R, créez un nouveau script R et préparez votre environnement de travail:

# Chargez la librairie `tidyverse` (aide: utilisez la fonction `library()`)


# Importer `burghardt_et_al_2015_expt1.txt` et mettez le dans un objet appelé `expt1`
#(aide: utilisez la fonction `read_tsv()`)


1. Ajouter des résultats de tests statistiques à des graphiques

Pour pouvoir rajouter des résultats de tests statistiques à des graphiques, nous allons utiliser des fonctions de la librairie ggpubr.

library(ggpubr)

Camparaison de moyennes de deux (ou plusieurs) groupes

Afin de faire un test de comparaison de moyennes, nous utilisons la fonction stat_compare_means.

ggplot(expt1, aes(genotype, days.to.flower, colour = fluctuation)) +
  geom_boxplot() +
  stat_compare_means(label = "p.format", method="t.test")
## Warning: Removed 83 rows containing non-finite values (stat_boxplot).
## Warning: Removed 83 rows containing non-finite values (stat_compare_means).

Comme nous avons deux groupes à comparer (température fluctuante vs température constante), la fonction fait la comparaison de ces deux groupes automatiquement pour chaque génotype et donne la p-value associée.


Plutôt que d’ajouter les valeurs des p-values, il est aussi possible d’utiliser un code à base d’étoiles pour indiquer le niveau de significatifité du test:

ns: p > 0.05

*: p <= 0.05

**: p <= 0.01

***: p <= 0.001

****: p <= 0.0001

ggplot(expt1, aes(genotype, days.to.flower, colour = fluctuation)) +
  geom_boxplot() +
  stat_compare_means(label = "p.signif", method="t.test")
## Warning: Removed 83 rows containing non-finite values (stat_boxplot).
## Warning: Removed 83 rows containing non-finite values (stat_compare_means).


Par contre, s’il n’y a pas deux groupes clairement identifiés, tous les groupes vont être comparés, en utilisant par défault avec un test de Kruskal–Wallis:

ggplot(expt1, aes(genotype, days.to.flower)) +
  geom_boxplot() +
  stat_compare_means(label = "p.format")
## Warning: Removed 83 rows containing non-finite values (stat_boxplot).
## Warning: Removed 83 rows containing non-finite values (stat_compare_means).

Voir ce site pour en savoir plus sur les différentes options de tests de comparaison de moyennes (comparaison de deux groupes et de multiples groupes).


Tests de corrélation

Il est aussi possible d’ajouter le résultat d’un test de corrélation (spearman ou pearson) à un scatterplot en utilisant la fonction stat_cor():

ggplot(expt1, aes(blade.length.mm, rosette.leaf.num)) +
    geom_point()+
  stat_cor()
## Warning: Removed 333 rows containing non-finite values (stat_cor).
## Warning: Removed 333 rows containing missing values (geom_point).


Par défault, un test de pearson est effectué. Pour utiliser un test de spearman, il faut utiliser method="spearman"

ggplot(expt1, aes(blade.length.mm, rosette.leaf.num)) +
    geom_point()+
  stat_cor(method="spearman")
## Warning: Removed 333 rows containing non-finite values (stat_cor).
## Warning: Removed 333 rows containing missing values (geom_point).


Si vous avez plusieurs groupes, stat_cor va calculer la corrélation pour chaque groupe.

ggplot(expt1, aes(blade.length.mm, rosette.leaf.num, colour = fluctuation)) +
    geom_point()+
  stat_cor()
## Warning: Removed 333 rows containing non-finite values (stat_cor).
## Warning: Removed 333 rows containing missing values (geom_point).

Il est aussi possible d’ajouter une ligne de regression (linéaire ou non), incluant un intervalle de confiance, au graphique avec geom_smooth(). Notez cependant, comme on peut le voir ci-dessous, que cette fonction fonctionne mieux pour des regressions lineaires et non complexes.

ggplot(expt1, aes(blade.length.mm, rosette.leaf.num, colour = fluctuation)) +
    geom_point()+
  stat_cor() +
  geom_smooth()
## Warning: Removed 333 rows containing non-finite values (stat_cor).
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 333 rows containing non-finite values (stat_smooth).
## Warning: Removed 333 rows containing missing values (geom_point).

Exercice 1:

Ajouter le test statistique le plus approprié aux deux graphiques suivant:

ggplot(expt1, aes(fluctuation, rosette.leaf.num)) +
  geom_violin() 
## Warning: Removed 95 rows containing non-finite values (stat_ydensity).

ggplot(expt1, aes(days.to.bolt, days.to.flower)) +
    geom_point() 
## Warning: Removed 83 rows containing missing values (geom_point).

> > BONUS: Ajoutez une ligne de régression (pour le graphique où cela est possible). Déplacez le résultat du test statistique pour qu’il soit au meilleur endroit sur le graphique. Qu’est ce que vous pouvez conclure pour chacun de ces graphiques?


2. Ajouter des statistiques descriptives à des graphiques

Créer des chaînes de commandes avec les pipes %>%

Les “pipes” (%>%) permettent de faire une séquence d’opération sur des données, sans avoir besoin de créer des objets intermédiaires (ou de faire des commandes imbriquées très compliquées)

Grace au symbole %>% pipe, nous pouvons créer une chaîne de commandes. Pour cela nous devons d’abord faire une commande et ajouter %>% à la fin de la ligne qui va utiliser le résultat de cette commande comme input pour la fonction à la ligne suivante. Nous allons utiliser les pipes à partir de maintenant pour créer des chaines de commandes.

Extraire les statistiques descriptives avec group_by() et summarise()

Parfois nous voulons résumer nos données dans une table plus petite et en extraire des statistiques descriptives (moyenne, médiane, nombre d’observations …).

Ce type d’opération peut être fait avec la combinaison de deux fonctions: group_by() et summarise().

Notez que group_by() ne change pas le format de la table de données. Cette fonction liste des lignes qui doivent être groupées. Nous pouvons ensuite utiliser summarise() pour extraire des statistiques descriptives de chaque groupe.

Par exemple, nous pouvons extraire la moyenne pour le temps de floraison de chaque génotype:

group_by(expt1, genotype) %>% 
summarise(mean.days.to.flower = mean(days.to.flower, na.rm = TRUE))
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 10 x 2
##    genotype   mean.days.to.flower
##    <chr>                    <dbl>
##  1 Col Ama                   55.8
##  2 Col FRI                   70.2
##  3 fca-6                     70.4
##  4 flc-3 FRI                 53.1
##  5 flk-1                     77.6
##  6 fve-3                     81.2
##  7 ld-1                      89.8
##  8 Ler-1                     55.6
##  9 prmt5 FRI                 88.7
## 10 vin3-4 FRI                95.7

L’output contient deux colonnes:

  • genotype qui est la colonne qui a servi à grouper les données

  • mean.days.to.flower qui est la colonne crée par la fonction summarise

Il n’y a que 10 lignes dans cette table, une par génotype.

Il est possible de grouper les données par plus d’une variable.

Par exemple nous pouvons mesurer la moyenne, la médiane et l’écart type pour chaque génotype aux différentes températures:

group_by(expt1, genotype, temperature) %>% 
summarise(mean.days.flower = mean(days.to.flower, na.rm = TRUE),
          sd.days.flower = sd(days.to.flower, na.rm = TRUE),
          median.days.flower = median(days.to.flower, na.rm = TRUE))
## `summarise()` regrouping output by 'genotype' (override with `.groups` argument)
## # A tibble: 20 x 5
## # Groups:   genotype [10]
##    genotype   temperature mean.days.flower sd.days.flower median.days.flower
##    <chr>            <dbl>            <dbl>          <dbl>              <dbl>
##  1 Col Ama             12             66.6          25.5                62  
##  2 Col Ama             22             46.0          20.5                51  
##  3 Col FRI             12             81.1          26.2                80  
##  4 Col FRI             22             57.8          25.9                59.5
##  5 fca-6               12             76.2          20.9                72.5
##  6 fca-6               22             63.9          34.6                40  
##  7 flc-3 FRI           12             64.7          25.2                46  
##  8 flc-3 FRI           22             43            17.5                50  
##  9 flk-1               12             87.8          27.9                87  
## 10 flk-1               22             61.2          29.2                47  
## 11 fve-3               12             96.1          26.3                92  
## 12 fve-3               22             51.5           5.96               49  
## 13 ld-1                12             99.8          27.8               100  
## 14 ld-1                22             72.9          36.5                60  
## 15 Ler-1               12             67.4          22.7                64  
## 16 Ler-1               22             45.2          18.3                51.5
## 17 prmt5 FRI           12            103.           22.0               102  
## 18 prmt5 FRI           22             67.2          26.6                62  
## 19 vin3-4 FRI          12            111.           40.0                83  
## 20 vin3-4 FRI          22             72.5          29.7                69

Il y a maintenant 20 lignes dans la table, car chaque génotype apparaît deux fois (12 et 22 degrés)

Une autre information utile que nous pouvons extraire est le nombre d’observation pour chaque groupe. Pour cela nous devons utiliser la fonction n(), dans summarise() qui compte le nombre de ligne pour chaque groupe.

Par exemple, pour connaitre le nombre d’observations pour chaque génotype:

group_by(expt1, genotype) %>% 
summarise(n.obs = n())
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 10 x 2
##    genotype   n.obs
##    <chr>      <int>
##  1 Col Ama      135
##  2 Col FRI      128
##  3 fca-6         64
##  4 flc-3 FRI    136
##  5 flk-1         60
##  6 fve-3         60
##  7 ld-1          60
##  8 Ler-1         68
##  9 prmt5 FRI    125
## 10 vin3-4 FRI   121

Attention: Quand vous utilisez la fonction group_by(), les lignes du tableau restent groupées en fonction de la variable utilisée. Les opérations suivantes vont utiliser ces groupes, ce qui peut poser problème. Pensez à utiliser la fonction ungroup() pour enlever les groupes quand vous avez fini avec group_by() et summarise()

Exercice 2:

Calculez la médiane et l’écart-type de blade.length.mm et total.leaf.length.mm pour chaque genotype aux différentes day.length. Calculez aussi le nombre d’observations de chaque groupe.


Insérer les statistiques descriptives à un graphique

Une autre possibilité est d’ajouter les statistiques descriptives à un graphique contenant les données.

Pour cela, nous devons:

  • Utiliser group_by() et summarise() pour extraire les statistiques descriptives

  • Utiliser une fonction de la famille *_join() pour les combiner avec nos données

  • Nous pouvons maintenant faire un graphique contenant les données et les statistiques descriptives.

Par exemple, prenons ce boxplot:

  ggplot( expt1, aes(genotype, rosette.leaf.num)) +
  geom_boxplot()
## Warning: Removed 95 rows containing non-finite values (stat_boxplot).

Si nous voulons y ajouter le nombre d’observations pour chaque groupe, nous utilisons le script suivant:

group_by(expt1, genotype) %>% 
  summarise(n.obs=n()) %>% 
  mutate(n.obs=paste("n =",n.obs)) %>% 
  full_join(expt1, by="genotype") %>% 
  ggplot( aes(genotype, rosette.leaf.num)) +
  geom_boxplot() +
  geom_text(aes(label=n.obs, x=genotype, y=0))
## `summarise()` ungrouping output (override with `.groups` argument)
## Warning: Removed 95 rows containing non-finite values (stat_boxplot).


Exercice 3

Faites un voilin plot de total.leaf.length.mm pour chaque génotype et ajoutez la médiane pour chaque groupe (avec un point coloré) ainsi que le nombre d’observation de chaque groupe

BONUS: Modifiez le graphiqe que vous venez de créer pour ajouter les comparaison de moyenne de chaque génotype avec Col Ama (comme sur le graphique ci-dessous)


Avec cette même méthode il est aussi possible d’ajouter des informations pour des groupes formés à partir de deux variables.

Par exemple, si nous voulons faire un boxplot rosette.leaf.num pour chaque génotype en fonction de la température et y ajouter le nombre d’observations des différents groupes, nous utilisons le script suivant:

group_by(expt1, genotype, fluctuation) %>% 
  summarise(n.obs=n()) %>% 
  mutate(n.obs=paste("n =",n.obs)) %>% 
  full_join(expt1, by=c("genotype", "fluctuation")) %>% 
  ggplot( aes(genotype, rosette.leaf.num, fill=fluctuation)) +
  geom_boxplot() +
  geom_text(aes(label=n.obs, x=genotype, y=-2),position=position_dodge(0.8), angle=45)
## `summarise()` regrouping output by 'genotype' (override with `.groups` argument)
## Warning: Removed 95 rows containing non-finite values (stat_boxplot).

3. Selectionner des lignes avec la fonction filter

Avec la fonction filter(), nous pouvons garder toutes les lignes de notre table qui correspondent à des plantes qui ont subies une vernalisation.

Tout d’abord, nous devons connaitre les différentes valeurs de la colonne vernalization. Comme nous pouvons voir, il y a deux options: ‘NV’ et ‘V’.

unique(expt1$vernalization)
## [1] "NV" "V"

(note: $ permet de sélectionner une colonne en particulier de la table)

Comme nous voulons garder les plantes qui ont subies une vernalisation, nous devons filtrer les données pour garder les lignes pour lesquelles il y a “V” dans la colonne vernalization:

filter(expt1, vernalization == "V")
## # A tibble: 330 x 15
##    plant_nb genotype background temperature fluctuation day.length vernalization
##       <dbl> <chr>    <chr>            <dbl> <chr>            <dbl> <chr>        
##  1        1 Col Ama  Col                 12 Con                 16 V            
##  2        2 Col Ama  Col                 12 Con                 16 V            
##  3        3 Col Ama  Col                 12 Con                 16 V            
##  4        4 Col Ama  Col                 12 Con                 16 V            
##  5        5 Col Ama  Col                 12 Con                 16 V            
##  6        6 Col Ama  Col                 12 Con                 16 V            
##  7        7 Col Ama  Col                 12 Con                 16 V            
##  8        8 Col Ama  Col                 12 Con                 16 V            
##  9        1 Col Ama  Col                 12 Con                  8 V            
## 10        2 Col Ama  Col                 12 Con                  8 V            
## # … with 320 more rows, and 8 more variables: survival.bolt <chr>, bolt <chr>,
## #   days.to.bolt <dbl>, days.to.flower <dbl>, rosette.leaf.num <dbl>,
## #   cauline.leaf.num <dbl>, blade.length.mm <dbl>, total.leaf.length.mm <dbl>

Nous pouvons utiliser les opérateurs suivant pour définir les conditions pour filtrer les données:

Opérateur Condition de sélection Exemple
< inférieur à filter(expt1, days.to.bolt < 20)
<= inférieur ou égal à filter(expt1, days.to.bolt <= 20)
> supérieur à filter(expt1, days.to.bolt > 20)
>= supérieur ou égal à filter(expt1, days.to.bolt >= 20)
== égal à filter(expt1, days.to.bolt == 20)
!= différent de filter(expt1, days.to.bolt != 20)
%in% est contenu dans filter(expt1, genotype %in% c("Col FRI", "Ler-1"))

Il est aussi possible de combiner plusieurs conditions de sélection avec les opérateurs suivant:

Opérateur Signification Exemple
& ET filter(expt1, days.to.bolt == 20 & genotype == "Ler-1")
| OU filter(expt1, rosette.leaf.num < 8 | rosette.leaf.num > 100)

Nous pouvons aussi identifier les données manquantes (NA) avec la fonction is.na() ou sa négation (en utilisant !):

Opérateur Signification Exemple
is.na() données manquante filter(expt1, is.na(rosette.leaf.num))
!is.na() donnée non manquante filter(expt1, !is.na(rosette.leaf.num))

source de l’image

Par exemple, nous pouvons sélectionner les plantes qui ont été vernalisées ET qui ont poussées avec une température fluctuante:

filter(expt1, vernalization == "V" & fluctuation == "Var")

Il est aussi possible de sélectionner les plantes qui ont poussées avec 8h de jours OU qui fleurissent tardivement:

filter(expt1, day.length == "8" | days.to.bolt > 85)

Exercice 2: Filtrez les données pour garder les plantes selon les 3 cas de figures suivant (indépendants les uns des autres):

  1. Plantes qui ne sont pas du background Ler et qui ont été traitées avec une température fluctuante.
  2. Plantes qui ont fleuries (bolt) en moins de 57 jours et qui ont moins de 40 feuilles de rosette
  3. Plantes du génotype fca-6 pour qui le blade.length.mm n’est pas manquant

4. Représenter des données temporelles (gather)

Pour représenter des données temporelles avec une ligne par individu mesuré, nous allons utiliser un autre jeu de données qui contient des transcriptomes effectués sur des plantes d’Arabidopsis thaliana ayant été traités pour 15 minutes, 1 heure ou 4 heures avec une température de 17°C ou 27°C.

Tout d’abord, chargeons les données et gardons les dans l’object RNAseq:

RNAseq <- read_tsv("../data/data_expression_cortijo2017.txt")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   gene = col_character(),
##   Log2FoldToZero_17c_15min = col_double(),
##   TLog2FoldToZero_17c_1hr = col_double(),
##   Log2FoldToZero_17c_4hr = col_double(),
##   Log2FoldToZero_27c_15min = col_double(),
##   Log2FoldToZero_27c_1hr = col_double(),
##   Log2FoldToZero_27c_4hr = col_double(),
##   cluster = col_character()
## )

L’objectif est d’obtenir un graphique de ce type:

Retour à la page d’accueil